Modeling of the Electrostatic Interaction and Catalytic Activity of [NiFe] Hydrogenases on a Planar Electrode

Hydrogenases are a group of enzymes that have caught the interest of researchers in renewable energies, due to their ability to catalyze the redox reaction of hydrogen. The exploitation of hydrogenases in electrochemical devices requires their immobilization on the surface of suitable electrodes, such as graphite. The orientation of the enzyme on the electrode is important to ensure a good flux of electrons to the catalytic center, through an array of iron–sulfur clusters. Here we present a computational approach to determine the possible orientations of a [NiFe] hydrogenase (PDB 1e3d) on a planar electrode, as a function of pH, salinity, and electrode potential. The calculations are based on the solution of the linearized Poisson–Boltzmann equation, using the PyGBe software. The results reveal that electrostatic interactions do not truly immobilize the enzyme on the surface of the electrode, but there is instead a dynamic equilibrium between different orientations. Nonetheless, after averaging over all thermally accessible orientations, we find significant differences related to the solution’s salinity and pH, while the effect of the electrode potential is relatively weak. We also combine models for the protein adsoption–desorption equilibria and for the electron transfer between the proteins and the electrode to arrive at a prediction of the electrode’s activity as a function of the enzyme concentration.


■ INTRODUCTION
Hydrogenases are a family of enzymes that catalyze the redox reaction of hydrogen. During the oxidation of hydrogen, two electrons are removed from the hydrogen molecule, generating two protons. In the reduction pathway, two electrons are added to two protons, generating molecular hydrogen: 1−3

2H
2e H 2 + + Hydrogenases are found in bacteria, archaea, and some eukarya, living in a wide spectrum of environments: aerobic, anaerobic, 1 or extreme conditions like hydrothermal vents. 4 Their role is to oxidize H 2 to H + or reduce H + to H 2 , thus providing a reversible route for energy conversion depending on the cell's necessities. Hydrogenases coupled to other enzymes play a major role in the fermentation of biological substances to CH 4 and in microbial phosphorylation, where H 2 can serve as an energy source in place of NADH. In other metabolic routes, hydrogenases generate molecular hydrogen as a subproduct of reductive processes. 1 Depending on the organism, hydrogenases can be found floating free in the cytoplasm, associated with the cellular membrane, or in the periplasm of the cell as part of a reaction chain involving other proteins. In eukaryotic cells, they are often located in specialized compartments. 1,5 Hydrogenases can be categorized into three major groups, depending on the composition of the catalytic center. The catalytic center of [NiFe] hydrogenases ( Figure 1) consists of a heteronuclear core containing these two metals coordinated with cysteines, one −CN group, and two −CO groups. 6 A subclass of this group collects the so-called [NiFeSe] hydrogenases, which contains the same heteronuclear core but with the substitution of cysteine by selenocysteine within the catalytic center. The second class comprises [FeFe] hydrogenases, containing two iron atoms within the catalytic center, which is connected to a [4Fe4S] cluster. A third class comprises [Fe]-only hydrogenases. These enzymes have the peculiarity of having a single Fe atom within their catalytic center. 1,7 In recent years, hydrogenases have captured the interest of researchers due to their possible applications in the field of renewable energies. 3 Great efforts have been made in order to understand their structures and catalytic mechanisms, and how they could be applied to energy conversion. One of these approaches has been the incorporation of hydrogenases on the anode in fuel cells to produce electricity, or on the cathode of (photo)electrochemical cells to produce hydrogen, avoiding the use of rare, expensive, and poison-prone metals such as platinum. 8−10 Some researchers have addressed the use of immobilized hydrogenases directly on electrodes made of a diversity of materials, testing their performance in titanium, 9,11,12 gold, 4,13,14 silver, 15 or graphite 16 with different degrees of success. Other studies have explored the possibility of increasing the effective surface of the electrodes using porous carbon materials 16−18 or nanoparticles. 13 The diverse degrees of success in the experimental studies indicate that several factors, other than the choice of the enzyme and the substrate, could be involved in the rate of reaction. 3,5 Several authors have pointed out 19−21 that the orientation of the enzymes on the surface of the electrode has a direct effect on the reaction rate. Indeed, a whole review has been published recently on this topic. 22 As illustrated in Figure  1, the best orientations are believed to be those in which the external iron−sulfur cluster is as close as possible to the electrode. This favors the transfer of electrons from the electrode to the catalytic center, through the array of iron− sulfur clusters that works as a conducting wire. Note that the catalytic center is often protected from the exterior of the protein, as it resides within a pocket that prevents direct electron transfer to it. 1 For less favorable orientations, the electron transfer rate would determine the overall rate for the production or consumption of hydrogen, reducing the effectiveness of the enzyme.
Depending on whether the functionalization of the electrode is achieved by physical adsorption or covalent attachment, the protein orientation can be considered to be either reversible, and therefore satisfying an equilibrium condition, or irreversibly fixed at the moment of attachment. In turn, physical adsorption generally involves hydrogen bonding, van der Waals, hydrophobic, and electrostatic interactions. Electrostatic effects can be a dominant effect for the orientation of immobilized proteins, 23−26 even though the other factors may have a strong effect on the overall (orientationally averaged) protein−electrode interaction energy. 27 Electrostatic interactions between a protein and a conducting surface have sometimes been described in terms of the dipole moment of the isolated protein (i.e., without the electrode) that may be readily extracted from molecular dynamics simulations. 23 Within this scheme, the protein should preferentially absorb with its dipole orthogonal to the surface. However, the dipole provides only a very rough description of a complex charge distribution, which may fail at close range. The orientation may thus depend also on local effects, connected to the presence of specific functional groups or "patches" on the surface of the electrode or of the enzyme. 28 We shall return to this point below, in the discussion of our results.
Electrostatic effects in biomolecules and electrolyte solutions can be modeled by a range of methods, each with its strengths and weaknesses in terms of simplicity, generality, and computational cost. 29 One approach that balances these contrasting requirements is the one based on the Poisson− Boltzmann equation, which describes the distribution of mobile ions and the electrostatic potential around one or more charged objects, such as a protein, an electrode, or their combination. 30,31 Several software tools solve this equation numerically, including for example APBS, 32 Delphi, 33,34 MEAD, 35 MIBPB, 36 AFMPB, 37 and TABI. 38 Here we have selected PyGBe, 39−41 due to its ability to calculate efficiently the electrostatic interaction between multiple bodies at close range.
The present study uses computational simulations with PyGBe to calculate the electrostatic component of the interaction between a [NiFe] hydrogenase and an electrode, as a function of orientation, pH, salinity, and electric potential.
The key question we attempt to answer is whether it is possible to immobilize and control the orientation of the enzyme on the surface of an electrode by tuning these variables, without modifying the chemistry of the electrode (e.g., by a selfassembled monolayer of functional molecules or by covalent bonding of the protein to it) or of the protein (e.g., by introducing selected point mutations of the amino acids at its surface, away from the active site). We also offer a prediction of the activity of the functionalized electrodes as a function of the enzyme concentration, based on simple assumptions about the adsorption equilibria and the electron transfer between the proteins and the electrodes.

■ SYSTEMS AND COMPUTATIONAL METHODS
The Protein. The hydrogenase used in this study is a [NiFe] hydrogenase isolated from Desulfovibrio desulf uricans ATCC 27774. 42 This hydrogenase can be found in the periplasm of the cell associated with the periplasmic tetraheme cytochrome c 3 , taking care of the first step toward recycling the chemical energy liberated during the redox reaction of hydrogen back to the cytoplasm. It is a heterodimer of 89 kDa, with two subunits with molecular masses of 62 and 27 kDa. It contains an array of three iron−sulfur clusters: two [4Fe4S] clusters and one [3Fe4S] cluster, respectively, denominated distal (or external), proximal, and medial clusters. The Ni and Fe atoms within the catalytic core are coordinated by one −CN group, two −CO groups, and four cysteines. 42 This specific hydrogenase was chosen for our study mainly because of the good resolution in the published crystal structure (1.80 Å).
The crystal structure of the [NiFe] hydrogenase was obtained from the Protein Data Bank 43 as a PDB file (1e3d). This file contains the crystallized structure of the hydrogenase in the form of a tetramer, together with the solvent molecules and ions associated with the hydrogenase at the moment of crystallization. Thus, we preprocess the PDB file with the pdbeditor software, 43 in order to isolate the system to be used for the calculations. This includes only the A chain (small subunit, 27 kDa) and the B chain (big subunit, 62 kDa), without the water molecules and ions.
The Physical Model. The calculations of the interaction between the protein and the electrode are based on the linearized Poisson−Boltzmann equation. The enzyme is assumed to be rigid, in a conformation identical to that extracted from the crystal. It is immersed in an implicit solvent, consisting of water and mobile salt ions. 29−31 Fixed point charges are arranged inside the protein, at the positions of the atoms. The interior and exterior of the protein are defined by the solvent-excluded surface (SES). The relative permittivity (or dielectric constant) considered for the inside of the protein is ϵ 1 = 4, while the solvent region has the relative permittivity of water (ϵ 2 = 80). The electrode is modeled as a conductor with the geometry of a rectangular cuboid with dimensions 250 × 250 × 10 Å 3 . In a preliminary set of calculations, we checked that the calculated interaction energies do not change significantly upon further increases in the size of the electrode. The electrode behaves as a metallic conductor and does not have an associated permittivity, but its electrostatic potential has a fixed value at all points of a fine mesh representing its surface. The model is illustrated schematically in Figure 2.
The model described in Figure 2 leads to a system of partial differential equations that were solved numerically with the PyGBe program. 39−41 PyGBe�pronounced "pig-bee"�is based on a library of routines written in CUDA and exploiting the parallelism afforded by graphical processing units (GPUs), which can be driven by user-modifiable Python scripts. It was explicitly developed to calculate the electrostatic interaction between multiple bodies. It uses a boundary element approach to obtain the electrostatic potential ϕ(r), by solving the Poisson equation inside the protein (region Ω 1 ) and the linearized Poisson−Boltzmann equation in the surrounding solvent (region Ω 2 ): The summation in eq 1 runs over all the protein's atomic charges, while κ is the inverse of the Debye−Huckel screening length, which depends on the overall concentration of small ions dissolved in the aqueous medium: where c j and q j = ez j are the concentrations and charges of the ions, k B T is the thermal energy, and e is the elementary charge. The ionic strength I coincides with the salt concentration, for a monovalent salt such as NaCl (z j = ±1). For a physiological solution with I = 0.15 mol/L, one has κ = 0.125 Å −1 at room temperature. In addition, we have also considered a salt-free aqueous solution with I = 0.0 mol/L (neglecting the small contributions coming from H + and OH − ions). Equations 1 and 2 are coupled, because ϕ(r) and the electric displacement [−ϵ r ∇ϕ(r)] must be continuous at the interface between Ω 1 and Ω 2 (i.e., on the protein's SES). In addition, on the electrode's surface (Γ e ) we adopted a Dirichlet boundary condition, whereby the potential takes a constant value ϕ e : while the electric displacement at the water−electrode interface (more precisely, its component along the normal direction n) gives the local charge density (unit per area). Unlike the potential, which is constant throughout the electrode, this may depend on the position r: where σ(r) is the induced surface electric charge density.
In the calculations we included two ion-exclusion layers (denominated "Stern layers" within this study), respectively surrounding the protein and the electrode. Each layer corresponds to a region with a thickness of 2.0 Å, which has the dielectric constant of water (ϵ 2 = 80) but a local ion concentration equal to zero. The purpose of this layer is to prevent excessive accumulation of positive/negative ions in regions with very negative/positive potentials. It is essentially an empirical correction for the assumption inherent to the Poisson−Boltzmann model of electrolyte solutions, without any short-range correlations related to the size of the ions (assumed to be pointlike). The distance between the surface of the electrode and the van der Waals surface of the closest atom was kept at a constant value of 4.1 Å, in order to avoid overlap between their Stern layers. Below this distance, the continuum Region Ω 1 corresponds to the protein, with fixed point charges (black dots). It is bounded by a surface Γ 1 and a Stern layer. Region Ω 2 corresponds to the electrolyte solution. Surface Γ e corresponds to the boundary of the electrode, also surrounded by a Stern layer. hypothesis that underlies the Poisson−Boltzmann equation might be questionable. Since we did not optimize it, include other contributions to the energy, or allow changes in the protein structure, our calculations are likely to underestimate the strength of the protein−electrode interactions.
To better understand how different factors could affect the interaction of the hydrogenase with the electrode, a computational experiment was designed considering the following factors: solution pH (5, 6, 7, 8, and 9), electric potential of the electrode ϕ e (−0.05, 0.00, and 0.05 V), and salt concentration in the solution (I = 0.15 or 0.0 M). Note that the pH is limited to a range of values where protein denaturation is not expected to occur, while the electrostatic potential of the electrode is limited by the range of validity of the linearized Poisson− Boltzmann approach (|eϕ e | < k B T, where e is the elementary charge and k B T is the thermal energy). Note that ϕ e is an "electric potential" and its value is zero at infinity. It is not an "electrode potential" measured against some reference electrode. 44 Electrostatic Adsorption Energies. To calculate the electrostatic interaction between the protein and the electrode, it is first necessary to assign the proper charge and radius to the atoms, according to a specific force field, and define the surface of the protein. For this purpose, a series of "pqr" files containing the Cartesian coordinates of the atoms, their charges, and van der Waals radii were created using the pbd2pqr software. 45 One of the problems in assigning the atomic charges in hydrogenases is the presence of nonstandard atomic groups, with transition metal ions within the iron−sulfur clusters and in the catalytic centers of these enzymes. We determined their charges by quantum chemical calculations on the iron−sulfur clusters and the catalytic center extracted from the enzyme, where the dangling chemical bonds were saturated with CH 3 groups. These atom selections were used to calculate the CHELPG charges, 46 based on single-point unrestricted density function theory (DFT) calculations using the ORCA software. 47,48 A Gaussian basis set (def2-SVP) 49 was selected to perform the calculations with a tight self-consistent-field option, using the PBE0 50 hybrid density functional to compute the exchange−correlation energy. We performed calculations for the catalytic center with a total charge of 0 and spin multiplicities of 1, 3, 5, 7, and 9; for the external cluster, a total charge of 0 and multiplicities of 1, 3, 5, 7, and 9; for the medial cluster, a total charge of +2 and multiplicities of 2, 4, 6, 8, and 10; and for the proximal cluster, a total charge of −3 and multiplicities of 1, 3, 5, 7, and 9. In each case, we adopted the CHELPG charges corresponding to the spin multiplicity with the lowest energy. All these data are reported in the Supporting Information.
The Amber 46 force field file, used by pdb2pqr to create the final pqr files, was modified by adding the calculated CHELPG charges. The pdb2pqr program was then run using standard settings for the "propka" method. 45 This assigns the protonation state of the protein's ionizable groups, depending on the value of the solution pH. This protonation state was assumed to be fixed, independently of the distance and orientation of the protein on the electrode. In principle, this assumption could be relaxed, with an increase in calculation time. 24 The total charge on the protein and the modulus of its dipole moment are reported in Table 1, for each pH value. Note that the dipole moment of an object with nonzero charge depends on the choice of pole for its evaluation. Our values have been calculated with respect to the center of charge of the protein.
Simulations in PyGBe require that information about the protein structure is transferred to a mesh representing its surface. The mesh files for our calculations were created from the pqr files using Nanoshaper. 51 The following settings were adopted, seeking a balance between cost and numerical accuracy of the calculations: "grid scale" = 2.0, "smooth mesh" = true, "probe radius" = 1.4, and "keep water shaped cavities" = true. The probe radius of 1.4 Å is the standard value for a water molecule. This is considered as a sphere rolling around the protein, thus generating its SES. Example input files for Nanoshaper can be found in the Supporting Information.
Having defined the atomic charges and the meshes representing the protein and electrode surfaces, the PyGBe program was used to solve the linearized Poisson−Boltzmann equation. 39−41 The solution provided by PyGBe consists of the values of the electrostatic potential and the normal component of its gradient, for each grid point of each surface. These can then be used to obtain the electrostatic component of the free energy of the system. According to our assumption of reversible adsorption equilibrium, the free energy determines the probability that the enzyme adopts a specific orientation on the electrode.
In our study we considered the electrostatic component of the free energy of the system as the total energy of the hydrogenase−electrode system. This is calculated by PyGBe as the sum of Coulomb, solvation, and surface contributions, according to the equation The Coulomb energy is calculated from the Coulomb interactions of all point charges: where ε 1 is the dielectric constant within the protein, ϕ Coulomb (r j ) is the Coulomb potential at the position charge q j due to the other charges (q i ), and the double summation runs over all the protein's charged atoms. This can be a large but constant term, being independent of the protein's orientation and distance from the electrode. The solvation energy is the energy contribution from the protein's surroundings: solvent polarization, charged surfaces, and possibly other proteins. It is calculated as where ρ is the charge distribution and ϕ reac = ϕ − ϕ Coulomb is the electrostatic potential that arises due to the solvent's reaction, that contains the contribution of the polarization of the solute and of the solvent. Again, the summation runs over all the protein's atomic charges. Finally, the surface energy is where the integral is performed over the electrode's surface and Q e is the net charge on it. The interaction energy is calculated by subtracting the values of the energy of the isolated electrode and the isolated protein from the total energy: The values of the total energies for the protein and the electrode were calculated separately with PyGBe, using the same meshes of the combined calculation. A negative value of the interaction energy means the adsorbed state (protein@ electrode) is energetically more stable than the desorbed state. The (θ i , φ i ) arguments appearing on the left-hand side of eq 10 emphasize that this energy depends on the protein orientation, to be discussed under Protein Orientations and Probabilities.
Protein Orientations and Probabilities. Using PyGBe, we calculated the total energy for different orientations of the hydrogenase on the electrode. As shown in Figure 3, the orientation is defined by an inclination or tilt angle θ (0°≤ θ ≤ 180°) and an azimuthal angle φ (0°≤ φ < 180°). 41 In order to have an even and exhaustive sampling of the orientation, the tilt angle was incremented in even steps of dθ = 10°, while the azimuth was sampled using variable steps of dφ = 360°/max{1, 36 sin(θ)} (i.e., using only one point when θ = 0°and 180°and 36 points when θ = 90°). In this way, the differential solid angle associated with a specific combination of θ and φ has a roughly constant value: dΩ = sin(θ) dθ dφ. With these settings, the total number of sampled orientations was M = 390.
The total energy can be converted into probability associated with the orientation, since the orientations with the lowest energy should be those most likely of occur. According to the Boltzmann equilibrium distribution 52−54 where the normalization constant N is given by Note that the contributions representing different orientations can be summed evenly in eq 10, because the differential solid angles associated with them are identical. The interaction energies can be used in place of the total energies in eqs 11 and 12, leading to the same probabilities.
Adsorption Equilibria. The probabilities of eq 11 depend on the relative energies of the adsorbed states, but they are independent of the overall adsorption energy defined by eq 10. Adding or subtracting a constant value to the total energies leads to the same probabilities, because of the normalization in eq 12. However, a change in the overall adsorption energy will have an effect on the degree of coverage of the electrode by the proteins, for a given protein concentration in solution. This can be described using Langmuir's theory of adsorption. 52−54 Langmuir's theory assumes that the proteins do not interact with each other, neither in solution nor with the electrode. This is certainly an idealization, but it is consistent with our absorption calculations, which consider a single protein on the electrode.
Let χ be the overall coverage of the electrode, defined as the fraction of its area covered with proteins (0 ≤ χ ≤ 1). We also define χ(θ i , φ i ) as the fraction of the electrode area covered by the proteins with a specific orientation, out of M possibilities. These are related by Clearly, χ(θ i , φ i ) should be proportional to the probabilities of eq 11. We may thus write The Langmuir adsoption equation relates χ to the protein's osmotic pressure Π in solution: Here K is an equilibrium constant for the overall protein adsorption. χ is proportional to Π (and to the protein concentration in solution, if this behaves ideally) when Π ≪ K −1 , but it saturates at 1 when Π ≫ K −1 . The latter situation corresponds to the formation of a protein monolayer on the electrode surface. K summarizes the effects of all the individual equilibria, between hydrogenases in solution (H sol ) and hydrogenases on the surface (H surface ): Each has its own equilibrium constant K(θ i , φ i ), related as follows to the interaction free energies of eq 10: 54 The square brackets collect quantities depending on the vibrational motion of the adsorbed proteins on the surface (q vib is a vibrational partition function) and the osmotic pressure and the chemical potential of the proteins in solution, at some reference concentration (Π 0 , μ 0 ). These would be difficult if not impossible to calculate. Therefore, we simply assume them to be constant (independent of protein orientation). The 1/M prefactor in eq 17 accounts for the loss in rotational entropy, which occurs when the protein passes from the solution state to the adsorbed one, whereby it adopts a specific orientation. The overall equilibrium constant appearing in eq 14 can be obtained as the sum of the individual, orientation-dependent constants: Note the consistency among eqs 13 and 18. All these equilibrium constants depend on pH, salinity, and electrostatic potential at the electrode.
Electron Transfer Rates and Currents. As stated in the Introduction, the enzyme orientation on the electrode affects its catalytic efficiency, 2 because the active site may not be equally accessible to the reactants and products (H 2 and H + ) and electron transport to/from the electrode may be more or less difficult. Here we concentrate on the latter, since this is often the rate-limiting step for the redox reactions. 1,22 Modeling accurately the electron transport within a protein 55 and between the protein and the electrode is a complex task, which is beyond the scope of this article. Instead, we adopt a simple model which is based on the notion that electron transport in hydrogenases occurs by classical hopping or quantum mechanical tunneling�the distinction is not so important in the context of this paper�between the electrode and the external iron−sulfur cluster, and from there to the active site though the other clusters.
In general, the rate of electron transfer between two sites (the electrode and the external cluster, in this case) decays exponentially with the distance r: 52 k Ce r ET = (19) where β determines the rate of decay and C absorbs all other factors affecting the electron transfer rate (e.g., free-energy difference between the electron donor and acceptor, the reorganization energy, and temperature, according to the Marcus theory of electron transfer 52,55 ). In our case, each protein orientation is characterized by a different distance, so that r = r(θ i , φ i ). The distance from the iron atom coordinated to histidine 187 to the surface of the electrode was measured for this purpose. We have extracted the decay constant β = 0.45 Å −1 from a series of values calculated by Petrenko and Stein, for a model of a similar hydrogenase on a graphene platelet (see data and plots in the Supporting Information). 56 We point out that the value β = 1.4 Å −1 has been used in other publications (see, e.g., ref 57), but this is not expected to be universal and may depend on the protein's secondary structure features. 52 A numerical comparison of the results from these two values is also given in the Supporting Information.
The rate of hydrogen conversion by an individual enzyme with a specific orientation would be difficult, if not impossible, to measure. However, from a practical point of view, the most important quantity is the overall hydrogen conversion rate, or the overall electric current. As illustrated in Figure 4, the proteins adsorbed on a flat electrode may have different orientations, reflecting the probability distribution defined above. The total electric current can then be estimated by summing over all orientations, multiplying their probabilities by their respective electron transfer rates. The current should also the proportional to the area of the electrode (A) and to its degree of coverage (χ). We thus arrive at the expression where J 0 is a reference current density, independent of protein concentration and electrode area: In eq 21 we have assumed that a 0 , the area occupied by one adsorbed protein, is independent of its orientation. This is reasonable, considering the near-spherical shape of this hydrogenase.
The previous equations allow us to obtain the dependence of the electric current (or the hydrogen conversion rate) on the solution osmotic pressure (or the protein concentration in solution). Both the currents and the osmotic pressures can only be given in arbitrary units (AU), because of the uncertainties on the prefactors entering eqs 17 and 21. Despite this limitation, we will be able to calculate the dependence of the current on experimentally controllable variables such as the solution pH, salt concentration, and electrode potential. These prefactors will be discarded in the calculations that follow, but they may be used as adjustable parameters when fitting experimental data.
Analysis and Postprocessing. It is often desirable to analyze the results of a calculation in the form of an image, to understand better what happens when the protein comes close to the electrode. To this purpose, we wrote a Python script to generate a Visualization Toolkit (vtk) file that could be visualized using Paraview software. 58 A vtk file combines information on the meshes that represent surfaces of the protein and the electrode and the electrostatic potential calculated for every triangle of those meshes.

■ RESULTS AND DISCUSSION
Minimum-Energy Orientations. Different protein orientations (390 of them, overall) were tested for each combination of experimentally controllable variables: pH, electrostatic potential on the electrode ϕ e , and ionic strength I. The orientations with the lowest interaction energies are presented in Table 2, together with their probabilities and other salient properties. According to eqs 10 and 11, the probabilities are normalized by summing over all θ and φ angles (for given pH, ϕ e , and I). Figure 5 shows plots of the calculated electrostatic potential, for two representative minimum-energy orientations.
The interaction energies are negative in all cases, indicating that the adsorbed state is stable with respect to the one with the protein away from the electrode. On average, the ionic strength appears to be the most important variable determining the strength of absorption, measured by |ΔE min |. This is mostly lower than 10 kJ mol −1 in the presence of salt (with only one exception, for ϕ e = +0.05 V at pH 9), but it ranges between 13 and 53 kJ mol −1 in the salt-free cases. This is understandable, due to the screening effect of the dissolved ions on electrostatic interactions. Many experiments have corroborated the hypothesis that the amount of adsorbed proteins on charged surfaces decreases as the ionic strength increases. 27,59−63 An enhanced adsorption is in principle an advantage for electrode functionalization. We point out that it might not be feasible in practice to work under salt-free or salt-poor conditions, as they would reduce the electrical conductivity of the solution and possibly also the stability of the enzyme. However, the salt-free case is still interesting from a fundamental perspective, as it sets an upper limit to the strength and to the spatial range of the electrostatic interactions that can be expected in these systems.
Boubeta et al. 24 studied the interaction of lysozyme and other proteins with a negatively charged surface. They suggest that when electrostatic interactions are the main factor that determines the adsorption of a protein, two mechanisms can be considered to describe them. The first one is charge regulation (CR), whereby the protonation state of an amino acid depends not only on the pH of the solution but also on the chemical environment on the amino acid, the presence of charged groups, and the local dielectric properties. 64,65 Placing the protein near a highly charged surface will modify the pK a of the amino acids close to it. 66, 67 The second mechanism is charge patches (CP), whereby the protein tends to maximize The table gives also, for these orientations, the adsorption energy (ΔE min ), the probability (P min ), the distance of the external iron−sulfur cluster from the electrode (r min ), the closest amino acid to the electrode (AA min ), the net charge of the 10 closest amino acids (Q10 min ), and the overall charge on the electrode (Qe min ). The reference current densities (J 0 , eq 21) and the equilibrium constant for adsorption (K, eq 18) depend on all possible orientations.
The Journal of Physical Chemistry B pubs.acs.org/JPCB Article the number of oppositely charged amino acids next to the surface.
Concerning the role of the electrode's potential ϕ e , one might expect that the adsorption of the protein is favored when the sign of the potential is opposite that of the protein's net charge. However, a number of studies have shown that proteins can adsorb on modified surfaces even when the charges of the protein and the surface have the same sign, possibly due to a high local concentration of amino acids with the opposite charge (i.e., by a CP-type mechanism). [59][60][61]68,69 Indeed, we find that the electrostatic potential on the electrode affects the minimum-energy orientation of the protein, but it does not produce a large change in the overall absorption energy. In many cases, the changes are not easy to interpret. One situation where the changes appear to be systematic is at pH 8. Here the (θ min , φ min ) angles change when ϕ e changes from −0.05 to +0.05 V, indicating a reorientation of the protein due to the potential. However, |ΔE min | changes by a relatively small amount, from 5 to 7 kJ mol −1 . This is at first surprising, considering that the protein has a net change of −12.4e (see Table 1), and the overall charge on the electrode changes from −360e to +361e. However, one should also remember that the charge on the electrode is actually distributed, being roughly proportional to its size. Hence, a large electrode implies a large overall charge, but this does not automatically translate into a strong local interaction.
It is interesting to compare the total charge on the electrode when ϕ e = 0.00 V, as a function of pH and ionic strength. In this case, the total charge is zero when the protein is at infinite distance from the electrode. When the protein approaches the electrode and adsorbs on it, it induces a negative charge when it is positive (at acidic pH, see Table 1) and a positive charge when it is negative (at basic pH, see Table 1). This effect is strongest in the salt-free case, so that the induced charge on the electrode compensates almost exactly the total charge of the protein (about ±15 at pH 5 and 8), giving an electrically neutral electrode−protein complex.
Our calculations do not include the CR mechanism, because the protein charges are assumed to be independent of its orientation and distance from the electrode. However, we may consider the possible role of the CP mechanism. To do so, we have identified the amino acid closest to the surface for each orientation and calculated the total charge on the "patch" formed by the 10 closest amino acids. These data are also given in Table 2. One recurring orientation is that with θ = 80°and φ = 131°, in which the amino acid closest to the surface is Phe 354, which is electrically neutral and nonpolar. However, if we consider the whole patch, we find a value of −3e, at pH ≥6. This explains the recurrence of this particular orientation in Table 2 (7 instances out of 20). There is only one other case with a larger charge of the patch (+4e), at pH 5 in pure water. In this case the closest amino acid is a positively charged Lys 194, and the absorption energy is relatively large (−42 kJ/ mol). Interestingly, this is the closest amino acid also at pH 6, with a slightly different orientation. However, now the overall charge of the patch is only +1 (that of Lys itself), and consequently the value of the absorption energy is significantly smaller (−13 kJ/mol). Note that the presence of the Lys residue close to the electrode occurs also in other entries of Table 2, indicating that it may be crucial for achieving strong absorption on the electrode.
The probabilities of the minimum-energy orientations (P min in Table 2) are always relatively low, typically less than 0.1. This indicates that this orientation never dominates the distributions of absorption angles, but there is instead an equilibrium involving several other orientations that are only slightly higher in energy. There are only three entries in Table  2 where P min > 0.1. Two of them occur in the salt-free cases, where a large charge of the patch favors a very strong adsorption. The third case occurs in the saline solution at pH 5, with a positive electrode potential. The net charge of the patch is zero, but interestingly the amino acid closest to the surface is again a Lys.
With respect to the effect of the pH, many studies report that protein adsorption is maximum at pHs near the isoelectric point (IEP) of the protein. 60,70,71 Some of these observations could be ascribed to the role of nonelectrostatic forces (e.g., van der Waals and hydrophobic interactions), which may lead to protein adsorption even in the presence of unfavorable electrostatic interactions. 27,57,72 These nonelectrostatic forces could be included in an approximate way by a solventaccessible-surface-area (SASA) model. 73 In view of the nearspherical nature of our hydrogenase, the surface area is independent of orientation, unless the interaction is so strong as to deform and disrupt the proteins' three-dimensional structures (possible but unlikely under many practical  Thus, an overall enhancement or reduction of the protein− electrode interaction would leave the equilibrium distribution essentially unchanged. We also point out that the potential at the electrode is zero or relatively low in our calculations (±0.05 V). According to some reports, strong electric fields produced by applying higher voltages could change the preferential orientation of the protein and possibly even disrupt it, either by physical denaturation or by chemical reaction with radicals formed at the electrodes. 21,74−77 These phenomena are clearly crucial for the operation of enzymebased electrodes, but their investigation is beyond the possibilities of the present computational approach. Orientational Distributions. We now turn to discuss the consequences of the enzyme orientation on its catalytic efficiency. According to our previous discussion, this is expected to depend on the rate of electron transfer from the electrode to the external iron−sulfur cluster. Thus, the two orientations shown in Figure 5 can be expected to show different activities, since the external iron−sulfur cluster is far from the electrode on the left-hand side and very close to it on the right-hand side. The distances are reported in Table 2, under the r min column. The shortest distance of all is 16.28 Å, obtained in a salt solution at pH 9 with a negative potential (θ = 130°and φ = 235°).
However, specific orientations are not representative of the whole system, because they account for a relatively small fraction of all the adsorbed proteins. Those in Figure 5 have probabilities of P min = 0.198 and 0.156, respectively, while the third one discussed above has P min = 0.012 (see again Table 2). In all cases, more than 80% of the proteins would adopt an arrangement different from the minimum. There is a thermal and dynamical equilibrium between a multitude of orientations, each with a different catalytic activity. Figures 6 and 7 illustrate the absorption energies (panel a), the Boltzmann probabilities (panel b), and the electron transfer rates (panel c), as two-dimensional heatmaps as a function of the tilt (θ) and azimuthal (φ) angles. Note that Figures 6c and 7c are identical since, in the approximation of eq 19, the electron transfer rates depend only on the electrode−cluster distance. Figures 6d and 7d display the product of the probabilities and currents, as this is the quantity that determines the average contribution of each orientation to the current density, according to eq 21. Clearly, a high current density is obtained only when there is a good "overlap" between the Boltzmann probabilities (panel b) and the electron transfer rates (panel c).
The final results for the reference current densities (J 0 ) are reported in the last-but-one column of Table 2. Despite the large differences in the protein orientations at the minima, and in the electron transfer rates at these minima, the overall currents are remarkably similar. The highest and the lowest values of the current density, occurring in the salt-free solutions, differ by only 1 order of magnitude. If we consider only the salt-containing solutions, the differences are even lower.
Absorption Equilibria and Currents. As a final step, we consider the effect of protein concentration or, equivalently, its osmotic pressure (Π). According to the Langmuir model, this affects the electrode coverage (χ) through the overall equilibrium constant K (see eq 15). The value of K may be dominated by a few strongly adsorbed orientations, or it may result from several, roughly equivalent ones. Figure 8 illustrates the behavior of the total currents at zero electrode potential in salty and salt-free solutions, respectively (panels a and b). In the former, there is little difference between the results at different pHs, confirming our earlier conclusions. The solution at pH 5 seems to be marginally better, because of the combination of a slightly higher J 0 (determining the saturation value of the current at large protein concentrations) and a slightly higher K (determining the position of the inflection point in the Langmuir isotherms). There is a much greater dependence on pH in the salt-free solutions. The cases at pH 8 and 9 are characterized by strong absorption (large K), but the enzyme is essentially locked in an unfavorable orientation for electron transfer (small J 0 ). At pH 6 and 7 we have the opposite condition. By far, the best situation occurs again at pH 5, thanks to a favorable combination of large K and large J 0 . Figure 9 demonstrates the effect of modulating the electrode potential, within the range compatible with the applicability of the linearized Poisson−Boltzmann equation (−0.05 V ≤ ϕ e ≤ 0.05 V). The most visible changes occur at the two ends of the pH spectrum, where the protein has either a large and positive charge or a large and negative charge. A negatively charged protein (pH 9) has a much greater affinity for the electrode at positive potential. Interestingly, in order of decreasing affinity, the electrode at negative potential comes before the electrode at zero potential. This is counterintuitive and could not have been expected on the basis of the protein's charge and dipole. Considering the positively charged protein (pH 5), absorption is more favorable with a negative potential, but the sensitivity of the equilibrium constant to this variable is much lower in this case. There is however a greater sensitivity of the reference current density, which can be ascribed to a significant change in the distribution of the protein orientations. Note that the moduli of protein charge and dipole moment are similar in these two cases (see again Table 1). Again, we believe that this effect of the electrode potential could not have been easily predicted, on the basis of these descriptors of the protein charge distributions. As a final remark, we point out again that the Langmuir model neglects protein−protein interactions and does account for the possible formation of multilayer structures. These have been observed for a number of proteins, especially at high applied voltages (of the order of ±1 V), depending also on their degree of conformational rigidity. 78

■ CONCLUSIONS
We have presented the results of Poisson−Boltzmann calculations for the adsorption of a hydrogenase on a planar conducting electrode. One of the main goals was to determine the possibility to control the orientation of the protein, exploiting electrostatic interactions. These are not the only ones, but they are expected to dominate at least in the absence of a specific functionalization of the electrode (e.g., to produce a covalent attachment). According to our calculations, the adsorption of the enzyme on the electrode depends on the solution's pH and salinity and the electrostatic potential at the electrode. We have found that the latter does not affect much the adsorption energy of the protein, but its sign has a clear effect on the orientation. The strongest orientations are obtained in the presence of charge patches on the protein's surface, whose presence depends on the solution pH. The overall interaction is largest in salt-free solutions, where electrostatic interactions are not screened by the presence of counterions.
The total current, which measures the overall rate of the redox reactions at the active sites of the adsorbed enzymes, is maximum when there is a good match between the Boltzmann probabilities of the orientations, their electron transfer rates, and the strength of the protein absorption on the electrode surface. In general, we do not find situations dominated by a single protein orientation, but the overall current depends on a whole population of absorbed states. However, we have observed that a number minimum-energy states are characterized by a Lys residue close to the electrode's surface. It would be interesting to further pursue and check this observation, by studying other hydrogenases or by introducing further Lys residues on the protein's exterior, at favorable positions with respect to the external [3Fe4S] cluster.
Despite a number of approximations, we believe that the present computational approach has demonstrated its usefulness, first and foremost because of the physical insights can it can generate. We are now applying it to other hydrogenases, starting from the thermal-and oxygen-tolerant ones that are more interesting for renewable-energy applications. On a parallel line, we are interested in extending it by seeking the nonlinear solutions of the full Poisson−Boltzmann equation, including also charge regulation effects, larger values of the electrode potential, and hydrophobic effects by a SASA model. We point out that SASA models are widely used for protein solvation, but their applicability to electrode surfaces cannot be taken for granted. The configurations generated by these calculations would be useful also as a starting point for molecular dynamics simulations, accounting for more specific nonelectrostatic interactions and for the proteins' flexibility. Indeed, some flexibility is known to be essential for the functionality of an enzyme, for example, to assist the diffusion of the substrates to and away from the catalytic site. The molecular dynamics results could also be used the other way round, as an input to further Poisson−Boltzmann-based analyses of electrostatic effects and free energies 73 and to the evaluation of electron transfer rates. 56,79 Over a longer time scale, the enzymes' reorientational motion, which can be induced for example by a large switch in the electrode's voltage, 57,80,81 could also be treated by a coupling a nonlinear Poisson−Boltzmann solver with a microhydrodynamic description of the surrounding fluid. 82 Of course, such refinements of the model should be stimulated and cross-checked by experiments, including ones based on modern in situ operando spectroscopies that provide information about local interactions at the protein−electrode interface. 83,84 ■ ASSOCIATED CONTENT